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ABSTRACT 



We numerically investigate the dynamics of orbits in 3D circumbinary phase-space 
as a function of binary eccentricity and mass fraction. We find that inclined circumbi- 
nary orbits in the elliptically-restricted three-body problem display a nodal libration 
mechanism in the longitude of the ascending node and in the inclination to the plane 
of the binary. We (i) analyse and quantify the behaviour of these orbits with reference 
to analytical work performed by Farago & Laskar (2010) and (ii) investigate the sta- 



bility of these orbits over time. This work is the first dynamically aware analysis of the 
stability of circumbinary orbits across both binary mass fraction and binary eccentric- 
ity. This work also has implications for exoplanetary astronomy in the existence and 
determination of stable orbits around binary systems. 
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1 INTRODUCTION 



2 METHODS AND TERMINOLOGY 



The recent discovery of a circumbinary disk around the mi- 



croquasar SS433 (Blundell et al. 2008J along with evidence 
that such disks may be dynamically coupled to the accre- 
tion and outflow from s uch systems (jArtymowicz &; Lubow 
T996] |Regos et al.|[2QQ5] |Doolin fe Blundell||2009| |Perez M.' 
& Blundell 2010) has led us to investigate the behaviour of 



orbits encompassing binary systems. 

Additionally, in this era of exoplanetary astronomy we 
are finding new planets almost every week, but only a mi- 



nority those so far discovered are in binary systems (e.g. Lee 



|et al.|2009| |Beuermann et al.||2011| |Qian et al.|2011| . With 
an unknown but significant fraction of stars confined to bi- 
naries, and with new technologies and methods to detect 
circumbinary planets (e.g. |Schwarz et al.|2011 ), it is crucial 
that we understand the dynamics and stability of circumbi- 
nary orbits. 

In the following sections we investigate a circumbinary 
nodal libration with the aid of finely time-sampled numerical 
studies. We compare our results to analytic work performed 
by Farago &; Laskar ( 2010| and investigate the accuracy and 
limits of their model. 

With a second suite of numerical simulations we then 
determine the stability of the librating circumbinary orbits 
as a function of binary eccentricity and mass fraction. 



2.1 Orbital Elements 

A general orbit in 3D is described intuitively by the Keple- 
rian orbital elements (a, e, z, W, w, v) illustrated by Fig- 
ure [I] The eccentricity e, semi-major axis a and true 
anomaly v describe the motion of a body in its orbital plane, 
whilst the inclination i, longitude of the ascending node W 
and argument of perihelion w describe the orientation of the 
orbital plane with respect to some reference plane and direc- 
tion. The relationship between these orbital elements and a 
body's state vector (position and velocity) may be found in 
Green ( 1985 ) amongst others. 



In this coordinate space all trajectories are represented 
uniquely, with a Keplerian orbit having the property of con- 
serving all quantities but the true anomaly v(t) : which de- 
scribes the exact position of a body on its orbital path. 

In the circumbinary case we define the reference plane 
as the natural plane of the internal binary, and the reference 
direction as the vector along the binary's line of apseQ In 
our studies we consider the osculating orbital elements of 
test particles — that is to say their instantaneous orbital 



i.e. the binary system's semi-major axis 



© 0000 RAS 



2 Doolin & Blundell 




Table 1. Sampling of binary eccentricity and mass fraction 



Orbital Element 


min 


max 


A 


eccentricity 





0.6 


0.1 


mass fraction 


0.1 


0.5 


0.1 



Table 2. Sampling of circumbinary phase space where 

<2b = binary semi-major axis and Tj~, = binary orbital period. 




W longitude of the 
ascending node 



y 

Orbit projected onto celestial sphere 



Figure 1. A representation of the Keplerian orbital elements. 



elements about the barycentre of the binary system — as 
these quantities are no longer conserved. 



2.2 Numerical setup 

We have performed 3D numerical simulations of initially cir- 
cular (e = 0) circumbinar}j^] test particles in the elliptically 
restricted three-body problem. These massless test particles 
orbit synthetic binary systems of total mass 1M© and semi- 
major axis ah — 1AU. We scale all distances to the binary 
semi-major axis ah, and all times to the orbital period of 
the binary system Tb. 

We apply a customised adaptive step-size fourth and 
fifth-order Runge-Kutta integrator to integrate suites of test 
particles as a function of binary eccentricity eh and mass 
fraction ah — mi /(mi + 7712) where mi ^ m^. 

We track the orbital elements of each test particle about 
the centre of mass of the binary during integration, out- 
putting time-lapsed snapshots to a database. We take ad- 
vantage of the speed, organisation and SQL functionality of 
the database to study and accurately fit to the behaviours 
of orbits which we present in further sections. 

Each test particle is also monitored for instability. Un- 
stable orbits are identified and removed during integration 



Orbital Element 


min 


max 


A 


semi-major axis a 


1.5cib 


10a b 


0.5a,b 


inclination i 





71 


tt/20 


longitude of the ascending node W 





2tt 


tt/2 


true anomaly v 





2tt 


tt/2 


simulation length and snapshot At 




10 4 T b 


10T b 



'P-type' ( |Dvorak et al.|l989 > 



where a test particle is perturbed sufficiently from its initial 
orbit to approach either star, or if it evolves onto an un- 
bound trajectory (e > 1). This paper is not concerned with 
the ultimate fate of any test particle that experiences a close 
encounter with either stellar body. These test particles are 
rejected from the simulations and the precise scattering is 
not computed. 

Post-simulation stability criteria are applied to identify 
and reject test particles which do not quite reach escape 
velocity. 



3 THE NODAL LIBRATION 
3.1 Suite of simulations 

Our first suite of simulations was designed to be extensively 
time-sampled to expose the dynamics of circumbinary or- 
bits. The binary eccentricity and mass- fraction parameter 
space that we explore is laid out in Table ^whilst the initial 
phase-space sampling of test particles is laid out in Table [2] 
One additional spherical shell of test particles at a semi- 
major axis 50ab was also integrated for 5 x 10 6 Tb. 



3.2 Description 

As a typical test particle is integrated over the course of a 
simulation its semimajor axis a and eccentricity e remain 
constant. We discover a nodal libration mechanism in the 
inclination i and longitude of the ascending node W of all 
inclined orbits. 

This is best visualised in a polar (i cos W, i sin W) slice 
through parameter space, which may be termed a surface of 
section, as illustrated by Figure [2] Each line in Figure [2] is 
traced out by an individual test particle over the course of 
a simulation. Figure [2] reveals four distinct populations that 
we now examine in turn. 
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"'-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 

i*cos(W)/pi i*cos(W)/pi 




i*cos(W)/pi 

Figure 2. The (i cos W, i sin W) surface of section for circumbinary orbits as a function of binary orbital eccentricity e^. 
Green: prograde (i < tt/2). Precession is clockwise. 
Blue: retrograde (i > tt/2). Precession is anti-clockwise. 

Red: island of libration centred at (i = tt/2, W = tt/2). Precession is anti-clockwise. 
Purple: Island of libration centred at (i = tt/2, W = —tt/2). Precession is anti-clockwise. 



© 0000 RAS, MNRAS 000, 000-000 



4 Doolin & Blundell 



3.2.1 Prograde orbits (green) 

An inclination i = corresponds to a coplanar prograde cir- 
cumbinary orbit. Orbits of inclination ^ i < 7r/2, whilst 
not necessarily coplanar, we refer to as prograde. The pro- 
grade region of phase space therefore extends out from the 
centre (i — 0) of the surface of section (Figure [5| and en- 
compasses all orbits of a similar behaviour. 

Prograde orbits exhibit a precession in the longitude of 
the ascending node W . This evolution in W produces clock- 
wise paths around the surface of section shown in Figure [2] 



3.2.2 Retrograde orbits (blue) 

An inclination i = tv corresponds to a coplanar retrograde 
circumbinary orbit. Orbits of inclination tv/2 < i ^ tt, whilst 
not necessarily coplanar, we refer to as retrograde. The ret- 
rograde region of phase space therefore extends inwards from 
the outer limit (i — tv) of the surface of section (Figure [5| 
and encompasses all orbits of a similar behaviour. 

Retrograde orbits exhibit a precession in the longitude 
of the ascending node W. But counter to the prograde orbits 
the retrograde evolution in W produces anti-clockwise paths 
around the surface of section (Figure [2]). 

We expect that the precession in W observed in close- 
to-coplanar prograde and retrograde orbits is due to a cou- 
pling between the specific angular momentum of test parti- 
cles on inclined orbits and the z angular momentum of the 
binary. Such a coupling would exert a torque on the test 
particle producing a precession in the ascending node, akin 
to gyroscopic precession. 



3.2.3 Islands of libration (red & purple) 

An inclination i — tv/2 corresponds to a circumbinary or- 
bit which is exactly perpendicular to the binary plane. Fig- 
ure [2] shows two very clear libration islands centred on 
i = 7r/2, W = ±7r/2. A test particle on an orbit within a 
region of libration has its inclination i and ascending node 
W coupled to precess about the centre of libration. For both 
regions of libration this precession is ant i- clockwise. 



3.3 The geometry of the surface of section 

The geometry of the (i cos W, i sin W) surfaces of section 
shown in Figure [2] reveal a dependence on the internal binary 
eccentricity. More specifically, the extent of the two regions 
of libration can be seen to scale with binary eccentricity. A 
circular binary eb = exhibits no libration islands, whereas 
for a binary of eb = 0.6 the libration mechanism is becom- 
ing the dominant behaviour in phase space. We quantify the 
extent of the libration regions in § |3.5| 

Inspection of surfaces of section across values of binary 
mass fraction ah and radius a lead us to conclude that the 
geometry is both mass fraction and radius independent. The 
period of the precession of each test particle however does 



show a strong dependance on ah and a, which we explore in 
§|T7| 



3.3.1 Kozai cycles 

The geometry of the polar (i,W) surface of section (Fig- 
ure [5| with its islands of libration appears similar to that 
of the Kozai mechanism. Kozai (1962) showed analytically 
that an inclined circumstellai^ orbit may experience an os- 
cillating exchange between inclination and eccentricity, and 
also a libration in the argument of perihelion. But whilst the 
circumstellar (i,w) plane may share similarities with our cir- 
cumbinary (i,W) plane, we are dealing with a very different 
regime of the elliptically-restricted three-body problem. 



3.3.2 Symmetry 

The (i cos W, i sin W) surface of section (Figure [2]) is a 2D 
projection of the surface of a unit sphere, where i corre- 
sponds to the polar angle and W to the azimuthal angle. 
The central point of the surface of section (i = 0) corre- 
sponds to a co-planar prograde orbit (see Figure [I]) and at 
this point the longitude of the ascending node W is unde- 
fined. Equivalent ly one may regard this i — point as sit- 
uated at the north pole of a sphere, with the south pole at 
i — tv (a coplanar, retrograde orbit). Points of i — tv/2 iden- 
tify the equator. This type of spherical projection is known 
in cartographic circles as the azimuthal equidistant projec- 
tion. 

This unit sphere essentially defines the direction of 
the specific angular momentum vector h of a test particle. 
We note three particular planes of symmetry through this 
sphere. These may be specified by considering components 
of a circumbinary orbit's specific angular momentum h 
along axes (i) parallel to the line of apses of the binary, (ii) 
perpendicular to the line of apses of the binar}Q and (iii) 
perpendicular to the plane of the binary. We show these 
three planes of symmetry in Figure [3] via the colours: 



green: i — tv/2 
blue: We{0,Tv} 
red: W G {tt/2, 3tv/2} 



3.4 Previous work 



i z = 

^||apses 







Confirmation that the features of Figure [2] are not numerical 
artefacts comes reports of from similar behaviour in Verrier 
&; Evans| ( |2009| . In that work the authors report discovering 
a counter-play between the Kozai mechanism and a new 
circumbinary libration whilst modelling orbits within the 
double binary system HD98800. 



An orbit about one star of a binary system. 'S-type' jDvorak] 
et al.|1989| >. 

yet remaining in the binary orbital plane 
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Figure 3. The separatrix and symmetries of the (i cos W, i sin W) 
surface of section, h = specific angular momentum. 



Following on from Verrier & Evans (2009) is the ex- 



cellent analytic paper of Farago & Laskar (2010). In this 



article the authors consider the elliptically restricted three- 
body problem and take advantage of the assumption that, 
in the circumbinary case, the displacement of the third body 
r3 is greater than the relative separation of the binary r2i- 

Under this r% ^> 7*21 approximation Farago & Laska^ex- 
pand the three body Hamiltonian to second order in r2i/r^ 
and then average over the orbit of the binary and the third 
body to obtain a time-averaged quadrupolar Hamiltonian. 

This Hamiltonian promises to be very accurate in the 
regime r2i/r^ <C 1, as higher order terms in r^ilrz will tend 
to zero at a faster rate than those of second order. This 
model should not be so good for orbits closer to the binary 
system, where higher order terms will play a more substan- 
tial role. In the following sections we investigate the accu- 
racy of Farago & Laskar s model by testing the predictions 
that it makes, and the limits at which their quadrupolar 
approximation becomes insufficient. 



3.5 Separatrix and critical angle 

The separatrix (Figure [3] black) is the boundary in phase- 
space between different modes of behaviour. In our circumbi- 
nary surface of section the separatrix takes the form of a 
triple figure-of-eight, or two circles intersecting at (W — 0, 
i = 7r/2) and (W = 7r, i = 7r/2), dividing the regions of 
behaviour outlined in § |3.2| above. The separatrix which we 
plot in Figure [3] is actually the path of a rare test particle 
which, due to our stepping integrator, sampled more than 
one region of behaviour. 

In Section |3.3| we mentioned that the geometry of the 
surface of section is predominantly dependent on binary ec- 




0.4 0.6 

binary eccentricity 



Figure 4. Orbital behaviour in the W = tt/2 plane as a function 
of binary eccentricity. 
Green: prograde (i < tt/2) 
Blue: retrograde (i > tt/2) 

Red: island of libration centred at (i = tt/2, W - 



tt/2) 



centricity, and here we quantify this. We note that each point 
on the surface of section (Figure [5]) defines a unique path, 
and that every path intersects the vertical axis (W — ±7r/2). 
We define a critical angle z cr it as the inclination i at which 
the separatrix crosses the positive vertical axis (W — +vr/2) 
in the region ^ i ^ tv/2. 

Subsequent to the discussion of symmetry in § |3.3.2| it 
follows that the three other intersections of the separatrix 
with the vertical axis are reflections of the critical angle z cr it 
defined in the region (W = 7r/2, ^ i ^ tt/2). 



3.5.1 Measuring the critical angle 

We have run an additional suite of simulations to find and 
extract the critical angle as a function of binary eccentricity. 
These simulations were run with one shell of test particles at 
radius 25ab, longitude of the ascending node W — tt/2 and 
with a high resolution in inclination, intervals of Ai — 71/ '80, 
at various values of binary eccentricity. 

For each orbit sampled we visually identify the be- 
haviour (§ [3]) as a function of inclination and binary ec- 
centricity to find upper and lower boundaries on the sep- 
aratrix. The results are plotted in Figure ^] preserving the 
colour scheme of Figure [2] 

The green, red and blue points in Figure [4] represent the 
behaviour of test particles at these locations. We only plot 
the sampled points that lie either side of the separatrix. We 
extract the simplest polynomial fit to accurately describe 
the critical angle (0 ^ i cr it ^ tt/2) as a function of binary 
eccentricity eb with the constraint that the fit must pass 
between every pair of green-red points. This is given by: 
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where 



0.5 + aeh ± be\ ± ce^ ± < 



a = -0.7138 ±0.0023 
b = 0.1021 ±0.0030 
c = 0.5264 ± 0.0038 
d = -0.3942 ± 0.0047. 



(1) 



In all cases the angle of the separatrix in the region 
(tt/2 ^ Zcrit ^ 7r) is at inclination tt — z C rit, as it should be 
by symmetry arguments. 



3.5.2 Comparison with Farago & Laskar (2010) 



The time-averaged quadrupolar model of Fa rago &; La skar 
(2010) predicts the critical angle (their equation 2.34) and 



hence the location of the separatrix as 



^crit = arcsm 



1 



l±4eS 



(2) 



Since we measure the critical angle at a radius of 25ab 
from the binary we expect the quadrupolar approximation 
to hold, and indeed we find an excellent agreement between 
our measurements and Equation [2] The predicted location 
of the separatrix lies between every pair of our experimental 
limiting points. 

We plot our data, fit and the Farago & Laskar prediction 
together in Figure [5] There is an almost exact agreement be- 
tween model and data at this radius. The largest divergence 
is at eb — > 1, which is an unphysical limit as the binary itself 
becomes unbound. 



3.6 Constant of motion 



Verrier & Evans ( 2009 ) proposed an integral of motion (see 
their equation 3) for the libration islands to be the compo- 
nent of the specific angular momentum of an orbit along the 
line of apses of the internal binary, which may be expressed 



h\\ 



h sin i sin W. 



(3) 



This makes intuitive sense as the libration islands are 
centred at (i = tt/2, W = ±7r/2), i.e. the points at which 
a test particle's angular momentum is exactly parallel or 
antiparallel to the line of apses of the internal binary (see 
Figure [5]). So for small deviations from these central points 
we find that h\\ apses is conserved. 

Unfortunately this model breaks down as we move out 
from the centre of libration. In Figure [6] (upper panel) we 
show the Farago &; Laskar constant of motion (Eq[3} over 
the course of our integration for example test particles in 
proximity to the centre of the island of libration (i — 7r/2, 



0.4 0.6 

binary eccentricity 



Figure 5. Our experimental fit to the critical angle of the sepa- 
ratrix (Eq|lJ (black) and |Farago fc Laskar| s analytic expression 
(Eq§ (pink). 
Green: prograde (i < tt/2) 
Blue: retrograde (i > tt/2) 

Red: island of libration centred at (i = tt/2, W = tt/2) 



W = tt/2) from our simulation of binary eccentricity eb = 
0.6 and mass fraction ah = 0.5. We observe that this sug- 
gested constant of motion becomes insufficient for test par- 
ticles of lower inclination, which sample phase space away 
from the centre of the libration. 

A more promising integral of motion is given by 
the time-averaged quadrupolar model of |Farago &; L askar 
( 2010| ), their equation 2.20, which we translate to orbital 
elements as 



cos 2 i - 



sin 2 z(5 sin 2 W — 1) 



(4) 



We see that this equation reduces to the square of Equa- 
tion [3] for values of (i — » 7r/2, W — >• ±vr/2) close to the centre 
of libration. But this model correctly predicts the paths of 
orbits across the entirety of the surface of section and for val- 
ues of binary eccentricity eb . We plot this constant of motion 
in the lower panel of Figure [6] for the same test particles as 
above. 

We use the constant of motion predicted by the time- 
averaged quadrupolar model of iFarago &; Laskar (2010) as a 



test of the limits of the quadrupolar approximation. For each 
simulated test particle we have 10 3 snapshots over the course 
of integration. We therefore calculate an instantaneous /ifl 
for each test particle at each snapshot and ask the question: 
'how constant is the constant of motion? 7 

For example — in our shell of test particles at radius 
50ab we find that a typical test particle experiences variation 
in /ifl of rsj 1.4 x 10 -5 . This is measured by taking the 
standard deviation a of the instantaneous measurements of 
/ifl over the course of a simulation. Since /ifl is of order 
unity, this equates to a typical error of 0.0014% at these 
large radii. 
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Figure 6. A plot of the suggested constants of motion for example 
test particles in proximity to the centre of the island of libration 
(i = 7r/2, W = 7r/2) from our simulation of binary eccentricity 
eb = 0.6 and mass fraction = 0.5. Upper panel: Verrier & 
|Evans| ( |2009t (Eq[3), lower panel: |Farago Laskar] ( |2010| > (Eq]4|K 

We now investigate how well the quadrupolar approxi- 
mation holds as we consider orbits closer to the binary. To 
do so we create histograms of ct(/ifl) for each shell of test 
particles that we sample (see Table [2|. In Figure [7] we show 
histograms of a (/ifl) for radii 3, 4, 5 and 6<2b- We also colour 
the histograms according to the simulated binary eccentric- 
ity e b . 

As we consider orbits closer to the binary (a < 10ab), 
we see that the quadrupolar approximation begins to break 
down. At 6<2b most orbits are accurate to the quadrupolar 
model to better than 0.5%, but as we move in to 3a b some or- 
bits experience deviation of > 1%. We observe similar small 
perturbations in a test particle's semi-major axis a and its 
initially zero eccentricity e. 

We also note that orbits around binary systems of 
higher eccentricity (Figure [7] purple) deviate significantly 
more than those around circular binaries (Figure [7] light 
blue), which indicates that the higher order terms of the 
Farago fc Laskar] Hamiltonian have a greater dependence on 
binary eccentricity. 
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a/a_b = 3.0 




0.000 0.005 0.010 0.015 0.020 

std( h_FL ) 

a/a_b = 4.0 




0.000 0.005 0.010 0.015 0.020 

std( h_FL ) 

a/a_b = 5.0 




0.000 0.005 0.010 0.015 0.020 

std( h_FL ) 

a/a_b = 6.0 




0.000 0.005 0.010 0.015 0.020 

std( h_FL } 

Figure 7. Histograms of ct(/ifl) for radii 3, 4, 5 and 6ab, coloured 
according to binary eccentricity eb- 
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3.6.1 Significance of a (/ifl) 

In Figure [8] we plot four typical examples of the variation in 
/ifl over the course of integration. These are selected from 
a random sample of 1000 such plots of test particles from 
the 3ab radius bin (Figure [7] upper panel). 

Numerical noise from a stepping integrator would in 
principle accumulate over time. In Figure [8] we show that 
the variation in /ifl is present from the start and constant 
in magnitude over the duration of a simulation. Values of 
c(/ifl) are therefore not biased by an accumulation of nu- 
merical noise towards the end of the simulation. 

This absence of bias is additionally supported by the 
appearance of structure in Figure [8] most prominently in the 
lowest panel, which is further evidence of the physical effect 
of higher order terms than those present in the quadrupolar 



Hamiltonian of [Farago &; Laskar (2010) 



3.7 Period of precession 

From our finely time-sampled At = 10Tb data (Table [5| we 
extract a period P for the precession of each stable test par- 
ticle. These range from as little as 50Tb (or 5 x At) to greater 
than the simulation duration, but good fits are obtained for 
the vast majority. 

In Figure [9] we plot test particle traces on the 
(i cos W, i sin W) surface of section for our simulation of bi- 
nary eccentricity eb = 0.5 and binary mass fraction ah = 
0.5, as a function of test particle radius a /at,, colouring the 
traces by precession period. 



Figure 10. An example fit of P oc a n . This is for a prograde 
orbit of inclination i = 7r/4 about a binary of eccentricity eb = 
and mass fraction = 0.5. Plotted is precession period P (in 
units of binary orbital period T) vs test particle semi-major axis 
a (in units of binary semi-major axis a^). 



We observe that the period of precession P is correlated 
with a test particle's orbital radius and proximity to a sepa- 
ratrix. Orbits which appear to be missing from Figure [9] are 
unstable, and hence are not plotted. We explore the topic 
of stability in much greater detail in the second part of the 
paper (§|4j). 



3.7.1 P oc a n 

As discussed above and demonstrated by Figure [5] the closer 
a circumbinary test particle orbits to the binary system the 
shorter is its period of precession P. We find that a powerlaw 
provides a very good fit to P as a function of distance from 
the binary, and so we fit P oc a n to all orbits sampled, where 
n is a free parameter. An example fit is shown in Figure [l0| 
The data point from our simulation at radius 50ab provides 
a very good constraint. 

Our data show a tight clustering about n — 3.5, as 
plott ed in Figure [TT| This is not within the range given 
by |Verrier fe Evans| ( |2009| ) of n = 3.37 ± 0.06, b ut does 
agree with the time- averaged quadrupolar model of [Farago 
& Laskar which predicts n of exactly 3.5. 
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e 0.5 alpha 0.5 | a = 3.0 




e 0.5 alpha 0.5 | a = 3.5 



e 0.5 alpha 0.5 | a = 4.0 



e 0.5 alpha 0.5 | a = 4.5 






e 0.5 alpha 0.5 | a = 5.0 








'C05lW)/pi 

e 0.5 alpha 0.5 | a = 10.0 






Figure 9. The (i cos W, i sin W) surface of section for circumbinary orbits about a binary system of eccentricity = 0.5 and mass 
fraction = 0.5 at varius radii (scaled by binary semi-major axis a^), each track coloured by log period in binary orbital periods T^. 
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pro and retrograde orbits 
polar orbits 




130 3.35 3.40 3.45 3.50 3.55 3.60 3.65 3.70 



Figure 11. A histogram of fitted n in P oc a n showing two dis- 
tinct populations: the polar orbits and the pro/retrograde orbits. 



But on closer inspection of Figure we notice that 
P oc a n appears to consist of two populations. One of these 
populations is comprised of the prograde and retrograde or- 
bits, which cluster around n — 3.53, and the other corre- 
sponds to the polar orbits, which cluster around n = 3.47. 



3.7.2 An analytic expression for period of precession 



The time-averaged quadrupolar model of |Farago &: La skar 
(2010) makes the following prediction (their equation 2.32): 



1 



T b 3tt a h (l - a h ) \a h J y/( 
where h = /ifl as defined in Eq[4] 



/ a \V 2 F(k 2 )(l-e 2 ) 2 
W A /(l-eg)(/i + 4eg) 



(5) 



k 2 = 



be 2 



h 



and F(k 2 ) may be defined in terms of the complete elliptical 
integral of the first kind, K(k 2 ), as 




io 1 io 2 io 3 io 4 io 5 io 6 io 7 

P measured 



Figure 12. A plot of precession period P in units of binary orbital 
period Tb: measured vs predicted (Farago & Laskar 2010) for all 
particles in our simulation. 



group consists of test particles from our radius 50<2b simula- 
tions, hence the larger periods, whereas the lower left group 
represents every test particle from the main body of our 
simulations (1.5<2b ^ a ^ 10ab, summarised in Table |5j). 

The scatter in period at the upper end of the groups of 
Figure [12] is due to these periods being significantly greater 
than the duration of the simulation, and hence badly fit- 
ted. The scatter at lower periods (P < 10 3 Tb) is due to 
deviations of the experimental values from the quadrupo- 
lar modal at low radii. But within the well behaved boxed 
region of Figure |12| 

-PfL / -^measured — 1.008 ±0.029. 

The results of our numerical investigation into the dy- 
namics of circumbinary orbits in the elliptically restricted 
three body problem have correlated well with the analytic 
work of Farago & Laskar (2010). We have also explored the 



limits of their model due to the quadrupolar approximation. 



where 



K(k 2 ) where k 2 < 1 
K(k~ 2 ) where k 2 > 1 



K(k 2 



rn/2 

Jo 



Vl-k 2 i 



We calculate the expected period (Eq[5| for each test 
particle that we sample, and in Figure [l2| we plot a straight 
one-to-one comparison between the predicted and measured 
periods of precession. 

Figure [12] shows two populations — the upper right 



4 STABILITY OF CIRCUMBINARY ORBITS 

Figure [9] shows test particle traces on the (i cos W, i sin W) 
surface of section for our simulation of binary eccentricity 
eb = 0.5 and binary mass fraction ah = 0.5, as a function 
of test particle orbital radius a /ah from the centre of mass 
of the binary. Each radius bin is equally sampled but we 
only plot orbits which remained stable throughout the sim- 
ulation. It can be seen that at a radius of only 3ab very few 
of the initial test particles are actually stable, but at greater 
distances from the binary the surface of section fills in. 
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Table 3. Sampling of circumbinary phase space where 

<2b = binary semi-major axis and = binary orbital period. 



Orbital Element 


min 


max 


A 


semi-major axis a 


^ 1.5ab 


^ 5<2b 


0.05a b 


inclination i 





7T 


tt/80 


longt. of the asc. node W 


tt/2 


tt/2 




true anomaly v 





2tt 


tt/3 


sim. length and snapshot At 




5 x 10 4 T b 


200T b 



Figure [9] shows that in the plotted simulation (e b = 0.5, 
a b = 0.5) the closest stable orbits to the binary are those 
at the centres of the islands of libration. These orbits are 
stable at a radius of 3a b , with regions closer to the separatrix 
becoming stable out to 4a b . As we move further away from 
the binary the retrograde orbits begin to acquire stability 
for radii > 5a b , before the prograde orbits finally become 
stable at > 6a b . 

The time-averaged quadrupolar model of | Far ago fe] 



Laskar (2010) predicts the dynamics of circumbinary orbits 
but can make no attempt to discern whether these orbits 
are viable. In the following sections we reveal characteristics 
which are due to resonances between the binary and test 
particle orbital periods which the Farago & Laskar model 
cannot explain due to the time-averaging carried out in its 
derivation. 



4.1 Suite of simulations 

We ran a suite of simulations to investigate the stability 
of circumbinary orbits. As discussed in § |3.5| the dynamics 
of the circumbinary phase space are such that every orbit 
crosses the W — ±7r/2 axis and there exists a symmetry 
which reflects W — +7r/2 onto W — —tt/2. We may there- 
fore narrow down the region of phase space which we sample 
to only one value of the longitude of the ascending node, 

W +7T/2- 

We run these simulations to 50,000 binary orbital peri- 
ods, sampling phase space to a good density in orbital radius 
a and inclination i, as laid out in Table [3] We sample across 
binary eccentricity and mass fraction as in the above sec- 
tions (see Table [I]). 



4.2 A measure of stability 



As discussed in § |2.2[ each test particle is monitored for 
instability. Unstable orbits are identified and removed dur- 
ing integration where a test particle is perturbed sufficiently 
from its initial orbit to approach either star, or if it evolves 
onto an unbound trajectory (e > 1). Post-simulation stabil- 
ity criteria are applied to identify and reject test particles 
which do not quite reach escape velocity. 

Whereas in section |3] we were concerned with the orbits 




10000 20000 30000 4O000 50000 

time in units of binary orbital period T 

Figure 13. A histogram of the escape time Escape for particles 
which sample unstable orbits. 



which survived the simulation, here we are more interested 
in those which don't. In Figure [13] we plot a histogram of 
the escape times £ e sca P e at which particles are rejected from 
the simulation. The vast majority of unstable particles are 
caught at the start of the simulation — of the approximately 
half a million unstable test particles, over half of these are 
caught within the first 1000 binary orbital periods, with the 
distribution tailing off steeply even in log-space. 

Since we sample each orbit from multiple initial values 
of the true anomaly v (Table [T]) we measure an orbit's long- 
term stability by the fraction of initial test particles which 
survive the simulated duration of 5 x 10 4 Tb. 

In Figure [14] we show a density plot of this measure 
of stability across our entire parameter space. The major 
axes of this figure correspond to the simulation parameters 
of binary eccentricity and mass fraction, whilst the minor 
axes correspond to orbital radius and inclination. This is 
a density plot where each pixel corresponds to a sampled 
orbit, and the transparency to our measure of stability — a 
darker colour indicates a more stable orbit. We preserve the 
colour scheme of Figure [2] 

We draw the reader's attention to the following features 
of Figure [14] as revealed by our exquisitely detailed simula- 
tions: 

(i) First, and very broadly speaking, orbits are more sta- 
ble at lower binary eccentricity et>. 

(ii) With equal generality, retrograde orbits (blue) appear 
to be the most stable, followed by librating orbits (red), and 
finally prograde orbits (green). The difference in radius of 
the innermost stable orbit across inclination can be as large 
as 2db- 

(iii) There are vertical striations of instability, most no- 
ticeable in the higher eccentricity simulations eb ^ 0.5, and 
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Figure 14. Orbital stability plotted as a function radius a/a^ and inclination i on the W = tt/2 axis, across binary eccentricity — mass 
fraction parameter space. Colours: 
Green: prograde (i < tt/2) 
Blue: retrograde (i > tt/2) 

Red: island of libration centred at (i = tt/2, W = tt/2) 
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predominantly at inclinations of i ~ and i ~ 7r/2. We hy- 
pothesize that these regions of instability are due to orbital 
resonances between a test particle and the binary. 

(iv) We note very thin horizontal pinnacles of instability 
in non-librating orbits 0.2 ^ eb ^ 0.4. These pinnacles are 
located at inclinations i « ty/4 and i w 37r/4, and extend up 
to 3ab into otherwise stable phase space. 

(v) More central to the pinnacles discussed above are 
wider peninsulas of instability in the librating region. These 
appear symmetrically either side of i — tt/2 in the librating 
region for simulations of eb ^ 0.2 and converge upon each 
other as eb —> 0.6. 

(vi) The horizontal pinnacles and peninsulas are a func- 
tion of binary mass fraction at>. These features do not ap- 
pear in the a\y = 0.5 simulations, and are magnified towards 
increasingly extreme values of c^b- 



4.3 Previous work 

There exist numerous papers in the literature in which au- 
thors investigate long-term orbital stability within the copla- 
nar circular restricted three body problem, both numerically 
and analytically, and with emphasis on circumstellar and 
circumbinary orbits. 

With improvements in computation power authors have 
relaxed the circular constraint in the problem to investi- 
gate eccentric binary systems ( Dvorak et al.|[l989| Holman 
"fe Wiegert||1999| |Musielak et al.||2005|). But only recently 



have we had the computation power to relax the coplanar 
constraint on the problem to investigate inclined orbits. 

Pilat-Lohinge r~~et al.| ( 2003| performed three dimen- 
sional numerical experiments to determine inclined stability 
but they did not e xplore the complex l ibration structure of 
the phase spac^] Pila t-Lohinger et al.| also considered only 
prograde inclinations up to 50° and systems of equal mass 
fraction (ah = 0.5). 

This paper presents the first dynamic-aware analysis 
of the stability of inclined circumbinary orbits throughout 
binary mass fraction — eccentricity parameter space. 



4.4 Escape time 

In a companion figure to the stability plot of Figure [14 
show the escape time of each unstable orbit in Figure 
Here we use the same axes as Figure [14] — with major axes 
corresponding to the simulation parameters of binary eccen- 
tricity and mass fraction, and minor axes corresponding to 
orbital radius and inclination. This is a density plot where 



5 The inclined simulations of |Pilat-Lohinger et al.| contain only 
test particles of initial longitude of the ascending node W = 
(Pilat-Lohinge r et al.| private communication). As such they did 
not sample any librating orbits, which do not intersect the W = 
or W = 7T axes. 



each pixel corresponds to an orbit sampled, and the trans- 
parency to the inverse of escape time l/t eS ca P e — a darker 
colour indicates a longer surviving orbit. We again preserve 
the colour scheme of Figure [2] Since we consider multiple 
test particles per sampled orbit we take an average for the 
escape time. 

In Figure [15] we find matching features to those in Fig- 
ure [14] as described above in § |4.2| But here we can also ob- 
serve how long-lived the unstable orbits are. For example, in 
the low binary eccentricity simulations eb < 0.2 we find that 
orbits are either very quickly unstable, or definitely stable. 
But for the simulations of higher binary eccentricity eb and 
more imbalanced binary mass fraction ab, the most strik- 
ing features are long-lived. The test particles which sample 
phase space at these extreme points are almost stable, and 
survive for times of order the simulated duration ~ 50, 000Tb 
(Table [3}. 

Test particles on unstable orbits further inside the un- 
stable region (closer to the binary) are quickly accreted onto 
the binary stars and/or ejected from the system. It is within 
this region of instability that inflows and outflows are likely 
to be important. 



4.5 Further considerations 

The phase space underlying our study of stability within 
the three-body problem is arguably significantly chaotic. We 
believe that the features which we report are real and sig- 
nificant to the consideration of matter around binaries, but 
we advise caution before cranking up the phase-space reso- 
lution. We should consider the effect of perturbations away 
from this idealised model, such as the volume of the bodies 
considered and the feedback effect on the internal binary of 
a third body of non-negligible mass. 



5 CONCLUSIONS 

Our simulations show that inclined circumbinary orbits in 
the elliptically-restricted three-body problem demonstrate 
three distinct families of behaviour: close-to-coplanar pro- 
grade (i ~ 0) and retrograde (i ~ tt) orbits precess in the 
longitude of the ascending node, whilst close-to-polar orbits 
(i rsj 7r/2 and W ~ zLtt/2) have their longitude of the as- 
cending node and inclination coupled to precess about the 
centre of an island of libration. 

We have extracted the critical angle z cr it (eb) of the sepa- 
ratrix along the critical W = tt/2 axis between these regions 
of behaviour as a function of binary eccentricity. 

We have shown that the analytic time-averaged 
quadrupolar model of Farago & Laskar ( 2010| ) provides an 
excellent description of the behaviours of circumbinary or- 
bits at radii ^ 50ab . We have also shown that their model be- 
comes inaccurate to greater than 1% at orbital radii ^ 5ab, 
and especially in cases of high binary eccentricity. 
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Figure 15. Inverse escape time 1/tescape as a function of radius a/a^ and inclination i on the W = tt/2 axis, across binary eccentricity 
— mass fraction parameter space. Colours: 
Green: prograde (i < tt/2) 
Blue: retrograde (i > tt/2) 

Red: island of libration centred at (i = tt/2, W = tt/2) 
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With the first 3D dynamic-aware analysis of the sta- 
bility of circumbinary orbits we have discovered that these 
orbits are surprisingly stable throughout binary mass frac- 
tion — eccentricity parameter space. Our detailed simula- 
tions have put numerical limits on this stability and revealed 
complex structure in orbital radius — inclination space. 

Our work shows that circumbinary phase-space is rich 
and dynamic, full of remarkable and stable orbits which do 
not behave simply. We should not presume any given binary 
system to lack a circumbinary component unless otherwise 
demonstrated. Such a component may be a source of obscu- 
ration, emission, inflow or outflow. 
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